function out = RunDynareppPruning_4th(filenameModel,opt)

global out
% ************************************************************************
% This function does the following:
%  1) solves a model from DYNARE++ up to fourth order
%  2) saves the solution in the form 
%       y_t   = g(x_t-1,u_t,sig)
%       x_t   = h(x_t-1,u_t,sig)
%       u_t   = eps_t
%     where v_t = [x_t-1,u_t]
%  3) simulates the model using a pruning scheme.
%  3) simulates the model using a pruning scheme. 
% 
% INPUTS:
% - filenameModel: a string for the name of the mod-file. 
%   For instance: filenameModel = 'RBCmodel';
%
% INPUTS (optional)
% - opt.filenameDynarerepp: a string indicating where the dynare++ file is saved
%   For instance: filenameDynarepp = 'C:\Dynare++';
% - opt.orderApp: the approximation order of the model (must be <= 4)
% - opt.numSim: the length of the simulated time series from the model
% - opt.seedNum: a number to initialize the random number generator
% - opt.ySelect: string of control variables. We only report the solution 
%   and simulate for these variables.
% - opt.deleteAll: 1 for deleting all the Mat-files from Dynare++ else zero
% - AntitheticShocks: 1 for simulating with antithetic shocks, else 0
%
% IMPORTANT: Running this function requires linking to the following three functions
%  1) DynarePP_Unfold_Matrices_4th.m
%  2) StateSpace_form_LinearInov_4th.m
%  3) Simulate_Pruning_LinearInov_4th.m
%
% OUTPUTS
% - gSteadyState: the steady state values of y_t
%   Dimension: ny * 1
% - opt.gv: the first order derivatives for y_t with respect to v_t = [x_t-1;u_t]
%   Dimension: ny * nv
% - opt.gvv: the second order derivatives for y_t with respect to v_t = [x_t-1;u_t]
%   Dimension: ny * nv * nv
% - opt.gss: the second order derivatives for y_t with respect to sig (the perturbation parameters)
%   Dimension: ny * 1
% - opt.gvvv: the third order derivatives for y_t with respect to v_t = [x_t-1;u_t]
%   Dimension: ny * nv * nv * nv
% - opt.gssv: the third order derivatives for y_t with respect to (sig,sig,v_t)
%   Dimension: ny * nv
% - opt.gsss: the third order derivatives for y_t with respect to sig 
%   Dimension: ny * 1
% - opt.gvvvv: the fourth order derivatives for y_t with respect to v_t = [x_t-1;u_t]
%   Dimension: ny * nv * nv * nv * nv
% - opt.gssvv: the fourth order derivatives for y_t with respect to (sig,sig,v_t,v_t)
%   Dimension: ny * nv * nv
% - opt.gsssv: the third order derivatives for y_t with respect to (sig,sig,sig,v_t)
%   Dimension: ny * nv
% - opt.gssss: the fourth order derivatives for y_t with respect to sig 
%   Dimension: ny * 1
%
% - opt.hSteadyState: the steady state values of x_t
%   Dimension: nv * 1
% - opt.hv: the first order derivatives for x_t with respect to v_t = [x_t-1;u_t]
%   Dimension: nv * nv
% - out.eta: the square of the co-variance matrix to the innovations in the
%   extended state space system. 
%   Dimension: nv * nu
% - opt.hvv: the second order derivatives for x_t with respect to v_t = [x_t-1;u_t]
%   Dimension: nv * nv * nv
% - opt.hss: the second order derivatives for x_t with respect to sig (the perturbation parameters)
%   Dimension: nv * 1
% - opt.hvvv: the third order derivatives for x_t with respect to v_t = [x_t-1;u_t]
%   Dimension: nv * nv * nv * nv
% - opt.hssv: the third order derivatives for x_t with respect to (sig,sig,v_t)
%   Dimension: nv * nv
% - opt.hsss: the third order derivatives for x_t with respect to sig 
%   Dimension: nv * 1
% - opt.hvvvv: the fourth order derivatives for y_t with respect to v_t = [x_t-1;u_t]
%   Dimension: nv * nv * nv * nv * nv
% - opt.hssvv: the fourth order derivatives for y_t with respect to (sig,sig,v_t,v_t)
%   Dimension: nv * nv * nv
% - opt.hsssv: the third order derivatives for y_t with respect to (sig,sig,sig,v_t)
%   Dimension: nv * nv
% - opt.hssss: the fourth order derivatives for y_t with respect to sig 
%   Dimension: nv * 1
%
% - sigma: the square root of the covariance matrix for the structural innovations
%   Dimension: nu * nu
%
% - opt.label_y: string with the names of the variables in y_t
% - opt.label_v: string with the names of the variables in v_t
% - opt.label_u: string with the names of the variables in u_t
%
% - opt.ySim: simulated sample path for the control variables. 
%   Dimension: ny * numSim
% - opt.xSim: simulated sample path for the state variables 
%   Dimension: nv * numSim
%%
% ************************************************************************
%% Inputs
% non-optional input 
if nargin < 1
    error('please specify the first input argument');
end

% Default options
defopt.filenameDynarepp = which('dynare++.exe');
defopt.orderApp         = 4;
defopt.numSim           = 1e4;
defopt.seedNum          = 1;
defopt.ySelect          = 'all';
defopt.deleteAll        = 1;
defopt.AntitheticShocks = 0;

% Check if default options are overwritten by opts
if ~isfield(opt,'filenameDynarepp')
    opt.filenameDynarepp  = defopt.filenameDynarepp;
end
if ~isfield(opt,'orderApp')
    opt.orderApp = defopt.orderApp;
end
if ~isfield(opt,'numSim')
    opt.numSim = defopt.numSim;
end
if ~isfield(opt,'seedNum')
    opt.seedNum = defopt.seedNum;
end
if ~isfield(opt,'ySelect')
    opt.ySelect = defopt.ySelect;
end    
if ~isfield(opt,'deleteAll')
    opt.deleteAll = defopt.deleteAll;
end        
if ~isfield(opt,'AntitheticShocks')
    opt.AntitheticShocks= defopt.AntitheticShocks;
end     

%% Computing the approximated solution using Dynare++
% Identify the current directory
pathfile = [pwd '\'];

% Solving the model up to first order 
command  = [opt.filenameDynarepp,' --ord 1',' --sim 0 --per 0 --ss-tol 10e-12 --no-centralize --no-irfs ',pathfile,filenameModel,'.mod'];
[status, result] = dos(command);
if status
   error(result)
end
% Results are saved in the file filename_1.mat
movefile([pathfile,filenameModel,'.mat'],[pathfile,filenameModel,'_1','.mat']);
 
if opt.orderApp > 1
    % Solving the model up to second order 
    command  = [opt.filenameDynarepp,' --ord 2',' --sim 0 --per 0 --ss-tol 10e-12 --no-centralize --no-irfs ',pathfile,filenameModel,'.mod'];
    [status, result] = dos (command);
    if status
       error(result)
    end
    % Results are saved in the file filename_2.mat
    movefile([pathfile,filenameModel,'.mat'],[pathfile,filenameModel,'_2','.mat']);
end
if opt.orderApp > 2
    % Solving the model up to third order 
    command  = [opt.filenameDynarepp,' --ord 3',' --sim 0 --per 0 --ss-tol 10e-12 --no-centralize --no-irfs ',pathfile,filenameModel,'.mod'];
    [status, result] = dos (command);
    if status
       error(result)
    end
    % Results are saved in the file filename_3.mat
    movefile([pathfile,filenameModel,'.mat'],[pathfile,filenameModel,'_3','.mat']);
end
if opt.orderApp > 3
    % Solving the model up to fourth order
    command  = [opt.filenameDynarepp,' --ord 4',' --sim 0 --per 0 --ss-tol 10e-12 --no-centralize --no-irfs ',pathfile,filenameModel,'.mod'];
    [status, result] = dos (command);
    if status
       error(result)
    end
    % Results are saved in the file filename_4.mat
    movefile([pathfile,filenameModel,'.mat'],[pathfile,filenameModel,'_4','.mat']);
end

% We also delete the follwing files generated by Dynare++
delete([pathfile,filenameModel,'_f','.m']);  %filenameModel_f. 
delete([pathfile,filenameModel,'_ff','.m']); %filenameModel_ff.m
delete([pathfile,filenameModel,'.dump']);    %filenameModel.dump

%% Transformation of the approximated solution
% We load and unfold the output from Dynare
[f_ss,fv,fvv,fss,fvvv,fssv,fsss,fvvvv,fssvv,fsssv,fssss,sigma,nv,nu,nx,nz,dyn_vars,dyn_state_vars,dyn_shocks] ...
    = Dynarepp_Unfold_Matrices_4th(filenameModel,opt.orderApp);

% We set up the alternative state space representation
[gSteadyState_tmp,gv_tmp,gvv_tmp,gss_tmp,gvvv_tmp,gssv_tmp,gsss_tmp,gvvvv_tmp,gssvv_tmp,gsssv_tmp,gssss_tmp,...
    hSteadyState,hv,hvv,hss,hvvv,hssv,hsss,hvvvv,hssvv,hsssv,hssss,eta,ny_tmp,label_y_tmp,label_v,label_u] = ...
        StateSpace_form_LinearInov_4th(f_ss,fv,fvv,fss,fvvv,fssv,fsss,fvvvv,fssvv,fsssv,fssss,sigma,nv,nu,nx,nz,dyn_vars,dyn_state_vars,dyn_shocks);

% We only use the selected variables in g(v)
if strcmp(opt.ySelect,'all') ~= 1 
     % Ensure that y_select is a column vector of strings
     [ny1,ny2] = size(opt.ySelect);
     if ny2 > ny1
         opt.ySelect = opt.ySelect';
         ny          = ny2;
     else
         ny          = ny1;
     end
     gSteadyState    = zeros(ny,1);
     gv              = zeros(ny,nv);
     gvv             = zeros(ny,nv,nv);
     gss             = zeros(ny,1);
     gvvv            = zeros(ny,nv,nv,nv);
     gssv            = zeros(ny,nv);
     gsss            = zeros(ny,1);
     gvvvv           = zeros(ny,nv,nv,nv,nv);
     gssvv           = zeros(ny,nv,nv);
     gsssv           = zeros(ny,nv);
     gssss           = zeros(ny,1);
     label_y         = cell(ny,1);
     for i=1:ny
        for j=1:ny_tmp
            if strcmp(cellstr(opt.ySelect(i,:)),cellstr(label_y_tmp(j,:))) == 1
                label_y(i,:)     = label_y_tmp(j,:);
                gSteadyState(i,:)= gSteadyState_tmp(j,:);
                gv(i,:)          = gv_tmp(j,:);
                gvv(i,:,:)       = gvv_tmp(j,:,:);
                gss(i,:)         = gss_tmp(j,:);                
                gvvv(i,:,:,:)    = gvvv_tmp(j,:,:,:);                
                gssv(i,:)        = gssv_tmp(j,:);                
                gsss(i,:)        = gsss_tmp(j,:);
                gvvvv(i,:,:,:,:) = gvvvv_tmp(j,:,:,:,:);                
                gssvv(i,:,:)     = gssvv_tmp(j,:,:);                
                gsssv(i,:)       = gsssv_tmp(j,:);
                gssss(i,:)       = gssss_tmp(j,:);
            end
        end
    end
else
    ny           = ny_tmp;
    label_y      = label_y_tmp;
    gSteadyState = gSteadyState_tmp;
    gv           = gv_tmp;
    gvv          = gvv_tmp;
    gss          = gss_tmp;
    gvvv         = gvvv_tmp;
    gssv         = gssv_tmp;
    gsss         = gsss_tmp;
    gvvvv        = gvvvv_tmp;
    gssvv        = gssvv_tmp;
    gsssv        = gsssv_tmp;
    gssss        = gssss_tmp;
end
clear gSteadyState_tmp gv_tmp gvv_tmp gss_tmp gvvv_tmp gssv_tmp gsss_tmp ...
     gvvvv_tmp gssvv_tmp gsssv_tmp gssss_tmp ny_tmp label_y_tmp

% Simulation of the model using the pruning method
% [ySim,vSim] = Simulate_Pruning_LinearInov_4th(gSteadyState,gv,gvv,gss,gvvv,gssv,gsss,gvvvv,gssvv,gsssv,gssss,...
%     hSteadyState,hv,hvv,hss,hvvv,hssv,hsss,hvvvv,hssvv,hsssv,hssss,eta,opt.numSim,opt.seedNum,opt.orderApp,opt.AntitheticShocks);
% ySim = ySim + repmat(gSteadyState,1,opt.numSim);
% vSim = vSim + repmat(hSteadyState,1,opt.numSim);
% 
% % We delete the Mat.files created by Dynare
% if opt.deleteAll == 1
%     delete([pathfile,filenameModel,'_1.mat']);        %filenameModel_1.mat
%     if opt.orderApp > 1
%         delete([pathfile,filenameModel,'_2.mat']);    %filenameModel_2.mat
%     end
%     if opt.orderApp > 2
%         delete([pathfile,filenameModel,'_3.mat']);    %filenameModel_3.mat
%     end
%     if opt.orderApp > 3
%         delete([pathfile,filenameModel,'_4.mat']);    %filenameModel_4.mat
%     end
% end

%% The output is packaged into the struct "out"
% The approximated solution
out.gSteadyState = gSteadyState;
out.gv           = gv;
out.gvv          = gvv;
out.gss          = gss;
out.gvvv         = gvvv;
out.gssv         = gssv;
out.gsss         = gsss;
out.gvvvv        = gvvvv;
out.gssvv        = gssvv;
out.gsssv        = gsssv;
out.gssss        = gssss;
out.hSteadyState = hSteadyState;
out.hv           = hv;
out.eta          = eta;
out.hvv          = hvv;
out.hss          = hss;
out.hvvv         = hvvv;
out.hssv         = hssv;
out.hsss         = hsss;
out.hvvvv        = hvvvv;
out.hssvv        = hssvv;
out.hsssv        = hsssv;
out.hssss        = hssss;
out.sigma        = sigma;
out.label_y      = label_y;
out.label_v      = label_v;
out.label_u      = label_u;
% The simulated time series
% out.ySim         = ySim;
% out.vSim         = vSim;
% The options used for the computations
out.opt          = opt;
end
